!********************************************************************
!  LNGS, 21.12.03
!  The reasonable results were obtained on 23.12.03
!********************************************************************
! Update for mirror protons
! LNGS 11.05.05
!********************************************************************
! Come back to mirror neutrons in November-December 2010
! LNGS, 02.12.10
!********************************************************************
! LNGS, 26.02.12 sum of two solutions with our or mirror sources being 0
!********************************************************************

Subroutine tst
   use ircoeff, only: ReadIRCoeff, GetbeeIR, GetPIR
   implicit none

   real(8):: Ep, E0, z, r
   real(8),dimension(2):: b
   real(8),dimension(2,6):: P

   integer:: ic, im

!========================================================================

   call ReadIRCoeff

   E0 = 1d18
   Ep = 1d19
   z=0.0d0
   r = dlog(Ep/E0)
   call GetbeeIR(r,z,b)
   write(*,'(A,2G12.5)') ' b=', b(1), b(2)

   call GetPIR(r,z,P)
   do ic = 1,2
      do im = 1,6
	     write(*,'(2I2,G12.5)') ic, im, P(ic,im)
	  enddo
   enddo
	     
   print*, ' Ciao!'
   stop

end Subroutine tst
!----------------------------------------------------------------

!----------------------------------------------------------------
Program Mirror2
   use msimsl , only: Dmolch, TIMDY
   use commirr, only: ReadParamMirr                &
                    , ido, npdesp, nxp, ldy, nstep &
                    , tmax, tgrid                  &
                    , tol, hinit                   &
                    , xbreak, dlgEeV, y, dlgE3J    &
                    , jbegin      

   use coeff, only: ReadPCMB 

   implicit none

   integer:: j, k, npdes, nx

   real(8):: t, tend, tbegin
   integer:: IHOUR, MINUTE, ISEC, NSEC, NSECN, DSEC

   external fcnbc, fcnut

!========================================================================

!--------------
!   call tst
!--------------

   call ReadParamMirr      ! read parameters of the problem
   call ReadPCMB           ! read coefficients

   npdes = npdesp
   nx    = nxp

   write(*,'(A,E10.3,A,I5)') ' tol=',tol, ' nxp=',nxp

!--------------------------------------------------------------------
! Start calculations 

   if(jbegin == 1) then
      y = 0d0
      call WriteFile(jbegin,tbegin)   
   else
      call ReadFile(jbegin)   
   endif

   CALL TIMDY (IHOUR, MINUTE, ISEC)
   NSEC = ISEC + 60*(MINUTE + 60*IHOUR)

   do j = jbegin, nstep-1
      t    = tgrid(j)
      tend = tgrid(j+1)
!=====
      print*, ' t=',t,' t_end=',tend
      call dmolch (ido, fcnut, fcnbc, npdes, t, tend, nx  &
         , xbreak, tol, hinit, y, ldy)
      call WriteFile(j+1,tend)
!=====
      CALL TIMDY (IHOUR, MINUTE, ISEC)
      NSECN = ISEC + 60*(MINUTE + 60*IHOUR)
      DSEC = NSECN - NSEC
      NSEC = NSECN
!*************  print results with step 50 *******************
      write (*,'(a,i3,a,f10.7,a,I8,a)') '  step=',j,' at t =', t, ' in ',DSEC,' sec'
      write (*,'(8e11.3)') (y(1,k),k=1,nx,50)
      write (*,'(8e11.3)') (y(2,k),k=1,nx,50)
!*****************  end print results  ***********************
   enddo ! j
    
   call dmolch (3, fcnut, fcnbc, npdes, t, tend, nx  &
         , xbreak, tol, hinit, y, ldy)
   write(*,*) ' The workspace is free'  
   write(*,*) ' Calculation is over!'

end Program Mirror2
!********************************************************************

!********************************************************************
! t stands for \tau = \ln(\frac(1+z_f)(1+z))
!-----------------------------------------
! Cour(1) - U_r''  I  Cmir(1) - U_m''
! Cour(2) - U_r'   I  Cmir(2) - Um'
! Cour(3) - U_r    I  Cmir(3) - U_m
! Cour(4) - Ppn_r  I  Cmir(4) - Ppn_m
! Cour(5) - Q_r    I  Cmir(5) - Q_m
!-----------------------------------------
! 
   Subroutine fcnut (npdes, x, t, u, ux, uxx, ut)
   use msimsl
   use commirr, only: z_f1, Osc
   use coeff  , only: AprCoefZ, Cour, Cmir
 
   implicit none

   integer, intent(in):: npdes
   real(8), intent(in):: t, x
   real(8), dimension(*):: u, ux, uxx, ut

   real(8):: z1, z13 
       
   z1 = z_f1*dexp(-t)
   z13 = z1*z1*z1

   call AprCoefZ(z1,x)

   if(Osc == 'O') then
      UT(1) = (Cour(3) - Cour(4))*U(1) + Cour(2)*UX(1) + Cour(1)*UXX(1) + Cour(5) + Cmir(4)*U(2)         
      UT(2) = (Cmir(3) - Cmir(4))*U(2) + Cmir(2)*UX(2) + Cmir(1)*UXX(2)           + Cour(4)*U(1) 
   elseif(Osc == 'M') then
      UT(1) = (Cour(3) - Cour(4))*U(1) + Cour(2)*UX(1) + Cour(1)*UXX(1)           + Cmir(4)*U(2)         
      UT(2) = (Cmir(3) - Cmir(4))*U(2) + Cmir(2)*UX(2) + Cmir(1)*UXX(2) + Cour(5) + Cour(4)*U(1) 
   elseif(Osc == 'N') then
      UT(1) = (Cour(3) - Cour(4))*U(1) + Cour(2)*UX(1) + Cour(1)*UXX(1) + Cour(5) + Cmir(4)*U(2)             
      UT(2) = (Cmir(3) - Cmir(4))*U(2) + Cmir(2)*UX(2) + Cmir(1)*UXX(2) + Cour(5) + Cour(4)*U(1) 
   else
      print*, ' Check where sources are switched on: O/M'
      stop 
   endif
   end Subroutine fcnut
!********************************************************************

!********************************************************************
!  Boundary condition

   Subroutine fcnbc (npdes, x, t, alpha, beta, gamp)
   use commirr, only: xmin, xmax
   implicit none

   integer:: npdes
   real(8):: x, t
   real(8), dimension(*):: alpha, beta, gamp

   if(x == xmin) then
      alpha(1) = 0d0    ! This says not to care about the
      beta(1)  = 0d0    ! left boundary conditon
      gamp(1)  = 0d0
      alpha(2) = 0d0    ! This says not to care about the
      beta(2)  = 0d0    ! left boundary conditon
      gamp(2)  = 0d0
   elseif(x == xmax) then
     alpha(1) = 1d0
     beta(1)  = 0d0
     gamp(1)  = 0d0
     alpha(2) = 1d0
     beta(2)  = 0d0
     gamp(2)  = 0d0
   else
      print*, ' Look in fcnbc'
   endif

   end Subroutine fcnbc
!********************************************************************

!********************************************************************
   Subroutine WriteFile(j,t)
   use commirr, only: nxp, ldy, dlgEeV, y, dlgE3J  &
              , chzf, chgam, chEmax, chevl, chXmir &
              , Osc 
   implicit none
   integer:: j, i, k
   real(8):: t
   character(1):: chj1
   character(2):: chj2
   character(20):: flnm

! 1d-29 is a factor for approximate normalization to HiRes
   do k=1,nxp
      do i = 1,ldy
         if(y(i,k).gt.0d0) then
            dlgE3J(i,k)=3d0*dlgEeV(k) + dlog10(y(i,k)) - 29.4d0
         else
            dlgE3J(i,k)=0d0
         endif 
      enddo
   enddo
!------------------------------------
   if( j.lt.10) then
      write(chj1,'(I1)') j
      chj2 = '0'//chj1
   else
      write(chj2,'(I2)') j
   endif

   flnm = 'Res\'//Osc//chj2//chzf//chgam//chEmax//chevl

!------------------------------------
   if( Osc .eq. 'N') then
      open(2,file=flnm//'.dat',status='unknown')
   else
      open(2,file=flnm//chXmir//'.dat',status='unknown')
   endif
      write(2,'(I3,1X,E15.6)') j, t
      write(2,'(3F12.6)') (dlgEeV(k),(dlgE3J(i,k),i=1,2),k=1,nxp)
   close(2)

   end Subroutine WriteFile
!********************************************************************

!********************************************************************
   Subroutine ReadFile(j)
   use commirr, only: nxp, ldy, y, dlgE3J          &
              , chzf, chgam, chEmax, chevl, chXmir &
              , Osc 
   implicit none
   integer:: j, i, k
   real(8):: foo
   real(8),dimension(nxp):: dlgEeV    ! 
   character(1):: chj1
   character(2):: chj2
   character(20):: flnm

!------------------------------------
   if( j.lt.10) then
      write(chj1,'(I1)') j
      chj2 = '0'//chj1
   else
      write(chj2,'(I2)') j
   endif

   flnm = 'Res\'//Osc//chj2//chzf//chgam//chEmax//chevl

!------------------------------------
   if( Osc .eq. 'N') then
      open(2,file=flnm//'.dat',status='unknown')
   else
      open(2,file=flnm//chXmir//'.dat',status='unknown')
   endif
      read(2,'(I3,1X,E15.6)') j, foo
      read(2,'(3F12.6)') (dlgEeV(k),(dlgE3J(i,k),i=1,2),k=1,nxp)
   close(2)

   do k=1,nxp
      do i = 1,ldy
         if( dlgE3J(i,k) <= 0d0) then
            y(i,k) = 0d0
         else
            y(i,k) =  1d1**(dlgE3J(i,k)-3d0*dlgEeV(k)+29)
         endif 
      enddo
   enddo

   end Subroutine ReadFile
!********************************************************************

